clear all
set maxvar 120000
use data, clear
xtset city_code year
global control1 "ln_pop ln_pergdp ln_gov indu road"
global absorb1 "a(city_code year)"
global vce "vce(cl city_code)"
*******************************************************************************
**********************Heteroscedasticity Test**********************************
*******************************************************************************
reg Frequency treat##post $control1
estat imtest, white 
reg Number treat##post $control1
estat imtest, white 
*******************************************************************************
*****************************Multicollinearity Test***************************
*******************************************************************************
reg Frequency treat post $control1
estat vif
reg Number treat post $control1
estat vif

*******************************************************************************
***********************Serial Correlation Test*********************************
*******************************************************************************
quietly xtreg Frequency did17 $control1 i.year, fe
predict e, resid
gen e_L1 = L.e
gen e_L2 = L2.e
reg e e_L1 e_L2, vce(cluster city_code)
test e_L2
drop e e_*

quietly xtreg Number did17 $control1 i.year, fe
predict e, resid
gen e_L1 = L.e
gen e_L2 = L2.e
reg e e_L1 e_L2, vce(cluster city_code)
test e_L2
drop e e_*


*******************************************************************************
**************************************stat**************************************
*******************************************************************************
logout, save(table2) excel replace: tabstat Frequency Number treat post $control1, stats(mean sd min p50 max n) c(s)


*******************************************************************************
*********************************base line**************************************
*******************************************************************************
cap erase table3.rtf
eststo clear
foreach i in Frequency Number{
	eststo: reghdfe `i' treat##post, noa $vce
	estadd local variable "no"
	estadd local city_fe "no"
	estadd local year_fe "no"

	eststo: reghdfe `i' treat##post, $absorb1 $vce
	estadd local variable "no"
	estadd local city_fe "yes"
	estadd local year_fe "yes"

	eststo: reghdfe `i' treat##post $control1, $absorb1 $vce

	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
}
esttab using table3.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe) keep(*1.treat* *1.post*)
eststo clear

*******************************************************************************
*********************************parallel trend********************************
*******************************************************************************
forvalue year = 2012/2023{
	g treat`year' = treat if year==`year'
	replace treat`year' = 0 if missing(treat`year')
}
foreach i in Frequency Number{
	reghdfe `i' treat2012 treat2013 treat2014 treat2015 treat2016 treat2017 treat2018 treat2019 treat2020 treat2021 treat2022 treat2023 $control1, $absorb1 $vce
	coefplot, keep(treat2012 treat2013 treat2014 treat2015 treat2016 treat2017 treat2018 treat2019 treat2020 treat2021 treat2022 treat2023) coeflabels(treat2012="2012" treat2013="2013" treat2014="2014" treat2015="2015" treat2016="2016" treat2017="2017" treat2018="2018" treat2019="2019" treat2020="2020" treat2021="2021" treat2022="2022" treat2023="2023") vertical yline(0) ytitle("") addplot(line @b @at) ciopts(recast(rcap)) scheme(s1mono) title("`i'") level(95) xlabel(, angle(-45)) xline(6,lp(dash_dot)) xtitle("Year") ytitle("Coefficient")
	graph save "Graph" "gph_`i'.gph", replace
}

graph combine gph_Frequency.gph gph_Number.gph,xsize(11) ysize(5) scheme(s1mono)altshrink iscale(1.1)
graph export "figure1.png", as(png) name("Graph") replace
erase gph_Frequency.gph
erase gph_Number.gph
 


*******************************************************************************
************************************placebo************************************
*******************************************************************************
use data, clear
xtset city_code year
keep city_code year post treat post Frequency Number $control1
drop if missing(treat)
sort city_code
save placebo, replace

use placebo, clear
keep city_code
duplicates drop
save placebo_treat, replace

forvalue i = 1/300{
	use placebo_treat, clear
	gen random_digit= runiform() //生成随机数
	sort random_digit  //按新生成的随机数排序
	g treat_rand = 1 if _n<=70
	replace treat_rand = 0 if missing(treat)
	drop random_digit
	sort city_code
	save random_treat

	use placebo, clear
	merge m:1 city_code using random_treat, nogen
	erase random_treat.dta
	global control1 "ln_pop ln_pergdp ln_gov indu road"
	global absorb1 "a(city_code year)"
	global vce "vce(cl city_code)"

	reghdfe Frequency treat_rand##post $control1, $absorb1 $vce
	g b_f = _b[1.treat_rand#1.post]
	g se_f = _se[1.treat_rand#1.post]
	g t_f = b_f/se_f
	g p_f = 2*ttail(e(df_r),abs(t_f))

	reghdfe Number treat_rand##post $control1, $absorb1 $vce
	g b_n = _b[1.treat_rand#1.post]
	g se_n = _se[1.treat_rand#1.post]
	g t_n = b_n/se_n
	g p_n = 2*ttail(e(df_r),abs(t_n))

	keep b_f se_f t_f p_f b_n se_n t_n p_n
	duplicates drop
	save placebo`i', replace
}
erase placebo.dta
erase placebo_treat.dta

clear all
forvalue i = 1/300{
	append using placebo`i'
	erase placebo`i'.dta
}
tw (scatter p_f b_f, yaxis(1) msize(vtiny)) (kdensity b_f, yaxis(2) kernel(gaussian)), scheme(s1mono) legend(region(fcolor(none) lpattern(blank))) legend(label(1 "P value") label(2 "Kernel density of Coefficient")) xtitle("β") ytitle("P value", axis(1)) ytitle("Density", axis(2)) title("Frequency") xline(0.133,lp(dash_dot)) xlabel(-0.2(0.1)0.2)
graph save "Graph" "gph_1.gph", replace
tw (scatter p_n b_n, yaxis(1) msize(vtiny)) (kdensity b_n, yaxis(2) kernel(gaussian)), scheme(s1mono) legend(region(fcolor(none) lpattern(blank))) legend(label(1 "P value") label(2 "Kernel density of Coefficient")) xtitle("β") ytitle("P value", axis(1)) ytitle("Density", axis(2)) title("Number") xline(0.121,lp(dash_dot)) xlabel(-0.2(0.1)0.2)
graph save "Graph" "gph_2.gph", replace
*graph combine gph_1.gph gph_2.gph,xsize(11) ysize(5) scheme(s1mono)altshrink iscale(1.1)
grc1leg2 gph_1.gph gph_2.gph, legendfrom(gph_1.gph) xsize(10) ysize(5) scheme(s1mono) iscale(1) altshrink
graph export "figure2.png", as(png) name("Graph") replace
erase gph_1.gph
erase gph_2.gph

clear all
use data, clear
xtset city_code year
global control1 "ln_pop ln_pergdp ln_gov indu road"
global absorb1 "a(city_code year)"
global vce "vce(cl city_code)"

cap erase table4.rtf
eststo clear
foreach i in All Solo{
eststo: reghdfe `i' treat##post , $absorb1 $vce
estadd local variable "no"
estadd local city_fe "yes"
estadd local year_fe "yes"

eststo: reghdfe `i' treat##post $control1, $absorb1 $vce
estadd local variable "yes"
estadd local city_fe "yes"
estadd local year_fe "yes"
}
esttab using table4.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe) keep(*treat* *post*)
eststo clear


*******************************************************************************
**************************************PSM**************************************
*******************************************************************************
clear all
use data, clear
xtset city_code year
global control1 "ln_pop ln_pergdp ln_gov indu road"
global absorb1 "a(city_code year)"
global vce "vce(cl city_code)"

cap erase table5.rtf
eststo clear

psmatch2 did17 $control1, out(Frequency) ate radius caliper(0.05)ties
pstest

eststo: reghdfe Frequency treat##post if _weight!=., $absorb1 $vce
estadd local variable "no"
estadd local city_fe "yes"
estadd local year_fe "yes"

eststo: reghdfe Frequency treat##post $control1 if _weight!=., $absorb1 $vce
estadd local variable "yes"
estadd local city_fe "yes"
estadd local year_fe "yes"

psmatch2 did17 $control1, out(Number) ate radius caliper(0.05)ties
pstest

eststo: reghdfe Number treat##post if _weight!=., $absorb1 $vce
estadd local variable "no"
estadd local city_fe "yes"
estadd local year_fe "yes"

eststo: reghdfe Number treat##post $control1 if _weight!=., $absorb1 $vce
estadd local variable "yes"
estadd local city_fe "yes"
estadd local year_fe "yes"

esttab using table5.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe) keep(*1.treat* *1.post*) title(Panel A: Caliper=0.05)
eststo clear

psmatch2 did17 $control1, out(Frequency) ate radius caliper(0.01)ties
pstest

eststo: reghdfe Frequency treat##post if _weight!=., $absorb1 $vce
estadd local variable "no"
estadd local city_fe "yes"
estadd local year_fe "yes"

eststo: reghdfe Frequency treat##post $control1 if _weight!=., $absorb1 $vce
estadd local variable "yes"
estadd local city_fe "yes"
estadd local year_fe "yes"

psmatch2 did17 $control1, out(Number) ate radius caliper(0.01)ties
pstest

eststo: reghdfe Number treat##post if _weight!=., $absorb1 $vce
estadd local variable "no"
estadd local city_fe "yes"
estadd local year_fe "yes"

eststo: reghdfe Number treat##post $control1 if _weight!=., $absorb1 $vce
estadd local variable "yes"
estadd local city_fe "yes"
estadd local year_fe "yes"

esttab using table5.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe) keep(*1.treat* *1.post*) title(Panel B: Caliper=0.01)
eststo clear

*******************************************************************************
***********************Spatially Adjacent *************************************
*******************************************************************************
cap erase table6.rtf
eststo clear
foreach i in Frequency Number{
	eststo: reghdfe `i' treat##post  if edge==1, $absorb1 $vce
	estadd local variable "no"
	estadd local city_fe "yes"
	estadd local year_fe "yes"

	eststo: reghdfe `i' treat##post $control1 if edge==1, $absorb1 $vce
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
}
esttab using table6.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe) keep(*1.treat* *1.post*)
eststo clear

*******************************************************************************
**************************************ddml*************************************
*******************************************************************************
clear all
use data, clear
xtset city_code year
global control1 "ln_pop ln_pergdp ln_gov indu road"
global control2 "c.ln_pop##c.ln_pop c.ln_pergdp##c.ln_pergdp c.ln_gov##c.ln_gov c.indu##c.indu c.road##c.road"
global vce "vce(cl city_code)"
global fe1 "i.city_code i.year"
global fe2 "i.city_code i.year##i.province_code"

cap erase table7.rtf
foreach i in Frequency Number{
	eststo clear
	set seed 42
	dis "$fe1"
	ddml init partial , kfolds(5)
	ddml E[Y|X]: pystacked `i' $fe1, type(reg) method(rf)
	ddml E[D|X]: pystacked did17 $fe1, type(reg) method(rf)
	ddml crossfit
	eststo: ddml estimate $vce
	estadd local variable "no"
	estadd local variable_2 "no"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
	estadd local province_year_fe "no"

	set seed 42
	dis "$control1 $fe1"
	ddml init partial , kfolds(5)
	ddml E[Y|X]: pystacked `i' $control1 $fe1, type(reg) method(rf)
	ddml E[D|X]: pystacked did17 $control1 $fe1, type(reg) method(rf)
	ddml crossfit
	eststo: ddml estimate $vce
	estadd local variable "yes"
	estadd local variable_2 "no"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
	estadd local province_year_fe "no"

	set seed 42
	dis "$control2 $fe1"
	ddml init partial , kfolds(5)
	ddml E[Y|X]: pystacked `i' $control2 $fe1, type(reg) method(rf)
	ddml E[D|X]: pystacked did17 $control2 $fe1, type(reg) method(rf)
	ddml crossfit
	eststo: ddml estimate $vce
	estadd local variable "yes"
	estadd local variable_2 "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
	estadd local province_year_fe "no"

	set seed 42
	dis "$fe2"
	ddml init partial , kfolds(5)
	ddml E[Y|X]: pystacked `i' $fe2, type(reg) method(rf)
	ddml E[D|X]: pystacked did17 $fe2, type(reg) method(rf)
	ddml crossfit
	eststo: ddml estimate $vce
	estadd local variable "no"
	estadd local variable_2 "no"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
	estadd local province_year_fe "yes"

	set seed 42
	dis "$control1 $fe2"
	ddml init partial , kfolds(5)
	ddml E[Y|X]: pystacked `i' $control1 $fe2, type(reg) method(rf)
	ddml E[D|X]: pystacked did17 $control1 $fe2, type(reg) method(rf)
	ddml crossfit
	eststo: ddml estimate $vce
	estadd local variable "yes"
	estadd local variable_2 "no"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
	estadd local province_year_fe "yes"

	set seed 42
	dis "$control2 $fe2"
	ddml init partial , kfolds(5)
	ddml E[Y|X]: pystacked `i' $control2 $fe2, type(reg) method(rf)
	ddml E[D|X]: pystacked did17 $control2 $fe2, type(reg) method(rf)
	ddml crossfit
	eststo: ddml estimate $vce
	estadd local variable "yes"
	estadd local variable_2 "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
	estadd local province_year_fe "yes"

	esttab using table7.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable variable_2 city_fe year_fe province_year_fe) keep(did17) title(`i')
eststo clear
}


*******************************************************************************
*********************************spatial effect*********************************
*******************************************************************************
use data, clear
xtset city_code year
global control1 "ln_pop ln_pergdp ln_gov indu road"
global absorb1 "a(city_code year)"
global vce "vce(cl city_code)"
foreach i in Frequency Number{
	drop if missing(`i')
}
xtset city_code year
keep Frequency Number $control1 city_code year lon lat city_code year treat post
duplicates drop
xtset city_code year
xtbalance, range(2012, 2023)
xtset city_code year

preserve
keep city_code lon lat
duplicates drop
g a = 1
reshape wide lon lat, i(a) j(city_code)
sort a
save a, replace
restore

preserve
keep city_code lon lat
duplicates drop
foreach i in city_code lon lat{
	rename `i' row`i'
}
g a = 1
merge m:1 a using a, nogen
erase a.dta
drop a
reshape long lon lat, i(rowcity_code) j(city_code)
geodist lat lon rowlat rowlon, generate(dis)
keep *code dis
reshape wide dis, i(city_code) j(rowcity_code)
drop city_code
save w, replace
restore

xtset city_code year
foreach i in Frequency Number $control1{
	by city_code: ipolate `i' year , gen(`i'_ip) epolate
	drop `i'
	rename `i'_ip `i'
}
xtset city_code year
spatwmat using w.dta,name(w) standardize

cap erase table8.rtf
eststo: reghdfe Frequency treat##post $control1, $absorb1 $vce
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
eststo: reghdfe Number treat##post $control1, $absorb1 $vce
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
eststo: xsmle Frequency treat##post $control1, wmat(w) model(sdm) type(both) fe noeffects nolog cl(city_code)
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
eststo: xsmle Number treat##post $control1, wmat(w) model(sdm) type(both) fe noeffects nolog cl(city_code)
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
esttab using table8.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe)
eststo clear

erase w.dta


*******************************************************************************
**************************Serial Correlation Test*****************************
*******************************************************************************
use data, clear
xtset city_code year
global control1 "ln_pop ln_pergdp ln_gov indu road"
global absorb1 "a(city_code year)"
global vce "vce(cl city_code)"
cap erase table9.rtf
eststo clear
foreach i in Frequency Number{
	eststo: reghdfe `i' treat##post l(1/2).`i', $absorb1 $vce
	estadd local variable "no"
	estadd local city_fe "yes"
	estadd local year_fe "yes"

	eststo: reghdfe `i' treat##post $control1 l(1/2).`i', $absorb1 $vce

	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
}
esttab using table9.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe) keep(*1.treat* *1.post*)
eststo clear


*******************************************************************************
**************************************mechanism*********************************
*******************************************************************************
clear all
use data, clear
xtset city_code year
global control1 "ln_pop ln_pergdp ln_gov indu road"
global absorb1 "a(city_code year)"
global vce "vce(cl city_code)"

cap erase table10.rtf
eststo clear
foreach i in agglomeration ploy number_of_center{
	eststo: reghdfe `i' treat##post $control1, $absorb1 $vce
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
}
esttab using table10.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe) keep(*1.treat* *1.post*)
eststo clear
cap erase table11.rtf
eststo clear
foreach i in tech tech_pop edu edu_pop{
	eststo: reghdfe `i' treat##post $control1, $absorb1 $vce
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
}
esttab using table11.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe) keep(*1.treat* *1.post*)
eststo clear
cap erase table12.rtf
eststo clear
foreach i in mean_zx max_zx{
	eststo: reghdfe `i' treat##post $control1, $absorb1 $vce
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
}
esttab using table12.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe) keep(*1.treat* *1.post*)
eststo clear

*******************************************************************************
**************************************Heterogeneity *********************************
*******************************************************************************

clear all
use data, clear
xtset city_code year
global control1 "ln_pop ln_pergdp ln_gov indu road"
global absorb1 "a(city_code year)"
global vce "vce(cl city_code)"

sum market_index, detail
g market_dummy = 1 if market_index>8.64036
replace market_dummy = 0 if market_index<=8.64036
replace market_dummy = . if missing(market_index)
cap erase table13.rtf
eststo clear
foreach i in Frequency Number{
	eststo: reghdfe `i' treat##post $control1 if market_dummy==1, $absorb1 $vce
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"

	eststo: reghdfe `i' treat##post $control1 if market_dummy==0, $absorb1 $vce
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
}
esttab using table13.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe) keep(*treat* *post*)
eststo clear
set seed 42
bdiff, group(market_dummy) model(reghdfe ln_per_co_f_in treat##post $control1, $absorb1 $vce) reps(300) detail
bdiff, group(market_dummy) model(reghdfe ln_per_co_n_in treat##post $control1, $absorb1 $vce) reps(300) detail

bys province_code: egen shenghui_code = min(city_code)
g shenghui = 1 if city_code==shenghui_code
replace shenghui = 0 if missing(shenghui)
cap erase table14.rtf
foreach i in Frequency Number{
	eststo: reghdfe `i' treat##post $control1 if shenghui==1, $absorb1 $vce
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"

	eststo: reghdfe `i' treat##post $control1 if shenghui==0, $absorb1 $vce
	estadd local variable "yes"
	estadd local city_fe "yes"
	estadd local year_fe "yes"
}
esttab using table14.rtf, se b(%9.3f) se(%9.3f) compress append ar2 star(* 0.10 ** 0.05 *** 0.01) nogap scalars(variable city_fe year_fe) keep(*treat* *post*)
eststo clear
set seed 42
bdiff, group(shenghui) model(reghdfe ln_per_co_f_in treat##post $control1, $absorb1 $vce) reps(300) detail
bdiff, group(shenghui) model(reghdfe ln_per_co_n_in treat##post $control1, $absorb1 $vce) reps(300) detail










 

